#------------------------------------------------------------------------------#
######################## Robustness Checks ###################################
#------------------------------------------------------------------------------#

rm(list = ls())

# Setting directory: adjust accordingly to 

setwd()

#------------------------------------------------------------------------------#
# Loading and installing necessary packages ##################################
#------------------------------------------------------------------------------#

# Installing required packages 

packages <- c('tidyverse', 'fixest', 'modelsummary', 
              'performance', 'car', 'psych', 'REdaS', 'knitr', 
              'lmtest', 'sandwich', 'boost')

# Checking if is installed (and install if not)

packages.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

rm(packages, packages.check)

#------------------------------------------------------------------------------#
# Loading data base  #########################################################
#------------------------------------------------------------------------------#

Base <- readRDS(file = "data/Base.rds")

#------------------------------------------------------------------------------#
# Renaming Variables  ########################################################
#------------------------------------------------------------------------------#

Base <- Base %>%
  mutate(
    bf = `Belief climate change`,
    cau = `Causes of climate change`,
    imp = `Impacts of climate change`,
    pi_lf = `Political Ideology (Left-Right)`,
    pi_cs = `Political Ideology (Conservatism-Progressive)`, 
    sk = `Subjective knowledge`,
    ok = `Objective knowledge`,
    sc = `Scientific consensus`,
    ts = `Trust in scientists`, 
    nep = `The New Ecological Paradigm (NEP)`,
    ii = `Individualism index`,
    ei = `Egalitarianism index`,
    pe = `Personal Experience (extreme weather events)`)

# Variable Dictionary
mydict <- c("pi_lf" = "Political ideology: Left",  
            "pi_cs" = "Political ideology: Progressive",
            "sk" = "Subjective knowledge", 
            "ok" = "Objective knowledge",  
            "sc" = "Scientific consensus",
            "ts" = "Trust in scientists",
            "nep" = "The New Ecological Paradigm (NEP)",
            "ii" = "Individualism worldview",
            "ei" = "Egalitarianism worldview",
            "pe" = "Personal experience (extreme weather events)",
            "FemaleYes" = "Female", 
            "Escolaridade_3High school or equivalent" = "Education: High school or equivalent", 
            "Escolaridade_3Undergraduate or more" = "Education: Undergraduate or more",
            "ReligionCatholic" = "Religion: Catholic", 
            "ReligionEvangelical Pentecostal or other evangelical" = "Religion: Evangelical Pentecostal or other evangelical", 
            "ReligionEvangelical Traditional" = "Religion: Evangelical Traditional",
            "ReligionOthers/No Relig." = "Religion: Others/No Relig.",
            "Income1 - 2 minimum wages" = "Income: 1 - 2 minimum wages", 
            "Income2 - 3 minimum wages" = "Income: 2 - 3 minimum wages",
            "Income3 - 5 minimum wages" = "Income: 3 - 5 minimum wages", 
            "Income5 - 10 minimum wages" = "Income: 5 - 10 minimum wages",
            "Income10 minimum wages or more" = "Income: 10 minimum wages or more", 
            "IncomeDo not know/ Prefer to not answer" = "Income: Do not know/ Prefer to not answer", 
            "age" = "Age (Years)", 
            "Color4" = "Race: Black",
            "(Intercept)" = "Const")

#------------------------------------------------------------------------------#
# 1. Multicollinearity  ############################
#------------------------------------------------------------------------------#

multicol <- lm(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
               Female + Escolaridade_3 + Religion + Income + age + Color4, 
             data = Base, 
             weights = ponde2)

multicol_table <- multicollinearity(multicol)

kable(multicol_table, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",7)))

#------------------------------------------------------------------------------#
# 2. Multiple Hypotheses Testing  ############################
#------------------------------------------------------------------------------#

existence <- lm(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
               Female + Escolaridade_3 + Religion + Income + age + Color4, 
             data = Base, 
             weights = ponde2)

causes <- lm(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
               Female + Escolaridade_3 + Religion + Income + age + Color4, 
             data = Base, 
             weights = ponde2)

impact <- lm(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
               Female + Escolaridade_3 + Religion + Income + age + Color4, 
             data = Base, 
             weights = ponde2)

index <- lm(index_std ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
               Female + Escolaridade_3 + Religion + Income + age + Color4, 
             data = Base, 
             weights = ponde2)


# F Test: complete model
summary(existence)
summary(causes)
summary(impact)
summary(index)

# F Test: Psychological variables
existence_psych <- linearHypothesis(existence, c("sk = 0", "ok = 0", 
                                           "sc = 0", "ts = 0", "nep = 0", "ii = 0",
                                           "ei = 0", "pe = 0"), 
                                 white.adjust = "hc1")

causes_psych <- linearHypothesis(causes, c("sk = 0", "ok = 0", 
                                           "sc = 0", "ts = 0", "nep = 0", "ii = 0",
                                           "ei = 0", "pe = 0"), 
                                 white.adjust = "hc1")

impact_psych <- linearHypothesis(impact, c("sk = 0", "ok = 0", 
                                           "sc = 0", "ts = 0", "nep = 0", "ii = 0",
                                           "ei = 0", "pe = 0"), 
                                 white.adjust = "hc1")

index_psych <- linearHypothesis(index, c("sk = 0", "ok = 0", 
                                         "sc = 0", "ts = 0", "nep = 0", "ii = 0",
                                         "ei = 0", "pe = 0"), 
                                white.adjust = "hc1")

f_psych <- bind_rows(existence_psych, causes_psych, impact_psych, index_psych) %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(Model = c("Existence", "Causes", "Impacts", "Belief Index")) %>%
  dplyr::select(Model, "Res.Df", "Df", "F", "Pr(>F)")

kable(f_psych, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",6)))

# F Test: Political variables

existence_pol <- linearHypothesis(existence, c("pi_lf = 0", "pi_cs = 0"), 
                 white.adjust = "hc1")

causes_pol <- linearHypothesis(causes, c("pi_lf = 0", "pi_cs = 0"), 
                 white.adjust = "hc1")

impact_pol <- linearHypothesis(impact, c("pi_lf = 0", "pi_cs = 0"), 
                 white.adjust = "hc1")

index_pol <- linearHypothesis(index, c("pi_lf = 0", "pi_cs = 0"), 
                 white.adjust = "hc1")

f_pol <- bind_rows(existence_pol, causes_pol, impact_pol, index_pol) %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(Model = c("Existence", "Causes", "Impacts", "Belief Index")) %>%
  dplyr::select(Model, "Res.Df", "Df", "F", "Pr(>F)")

kable(f_pol, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",6)))

# F Test: Socio-demographic variables

existence_soc <- linearHypothesis(existence, c("FemaleYes = 0", "Escolaridade_3High school or equivalent = 0",
                           "Escolaridade_3Undergraduate or more = 0", "ReligionCatholic = 0",
                           "ReligionEvangelical Pentecostal or other evangelical = 0",
                           "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                           "Income1 - 2 minimum wages = 0", "Income10 minimum wages or more = 0",
                           "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                           "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                           "age = 0", "Color4 = 0"), 
                 white.adjust = "hc1")

causes_soc <- linearHypothesis(causes, c("FemaleYes = 0", "Escolaridade_3High school or equivalent = 0",
                           "Escolaridade_3Undergraduate or more = 0", "ReligionCatholic = 0",
                           "ReligionEvangelical Pentecostal or other evangelical = 0",
                           "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                           "Income1 - 2 minimum wages = 0", "Income10 minimum wages or more = 0",
                           "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                           "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                           "age = 0", "Color4 = 0"), 
                 white.adjust = "hc1")

impact_soc <- linearHypothesis(impact, c("FemaleYes = 0", "Escolaridade_3High school or equivalent = 0",
                           "Escolaridade_3Undergraduate or more = 0", "ReligionCatholic = 0",
                           "ReligionEvangelical Pentecostal or other evangelical = 0",
                           "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                           "Income1 - 2 minimum wages = 0", "Income10 minimum wages or more = 0",
                           "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                           "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                           "age = 0", "Color4 = 0"), 
                 white.adjust = "hc1")

index_soc <- linearHypothesis(index, c("FemaleYes = 0", "Escolaridade_3High school or equivalent = 0",
                          "Escolaridade_3Undergraduate or more = 0", "ReligionCatholic = 0",
                          "ReligionEvangelical Pentecostal or other evangelical = 0",
                          "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                          "Income1 - 2 minimum wages = 0", "Income10 minimum wages or more = 0",
                          "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                          "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                          "age = 0", "Color4 = 0"), 
                 white.adjust = "hc1")

f_soc <- bind_rows(existence_soc, causes_soc, impact_soc, index_soc) %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(Model = c("Existence", "Causes", "Impacts", "Belief Index")) %>%
  dplyr::select(Model, "Res.Df", "Df", "F", "Pr(>F)")

kable(f_soc, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",6)))

# F Test: Insignificant variables

existence_null <- linearHypothesis(existence, c("ii = 0", "ts = 0", 
                                                "pi_lf = 0", 
                                                "FemaleYes = 0", "Escolaridade_3High school or equivalent = 0",
                                                "Escolaridade_3Undergraduate or more = 0", "ReligionCatholic = 0",
                                                "ReligionEvangelical Pentecostal or other evangelical = 0",
                                                "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                                                "Income10 minimum wages or more = 0",
                                                "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                                                "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                                                "age = 0"), 
                                    white.adjust = "hc1")

causes_null <- linearHypothesis(causes, c("sk = 0", "ei = 0", "pe = 0",
                                          "pi_lf = 0",
                                          "FemaleYes = 0", "ReligionCatholic = 0",
                                          "ReligionEvangelical Pentecostal or other evangelical = 0",
                                          "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                                          "Income10 minimum wages or more = 0",
                                          "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                                          "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                                          "age = 0", "Color4 = 0"), 
                                   white.adjust = "hc1") 

impact_null <- linearHypothesis(impact, c("sk = 0", "ts = 0", "ei = 0", "pe = 0", 
                                          "pi_cs = 0",
                                          "FemaleYes = 0", "ReligionCatholic = 0",
                                          "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                                          "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                                          "age = 0", "Color4 = 0"), 
                                white.adjust = "hc1")

index_null <- linearHypothesis(index, c("sk = 0","ei = 0", 
                                        "pi_cs = 0",
                                        "FemaleYes = 0", "ReligionCatholic = 0",
                                        "ReligionEvangelical Traditional = 0", "ReligionOthers/No Relig. = 0",
                                        "Income1 - 2 minimum wages = 0", "Income10 minimum wages or more = 0",
                                        "Income2 - 3 minimum wages = 0", "Income3 - 5 minimum wages = 0",
                                        "Income5 - 10 minimum wages = 0", "IncomeDo not know/ Prefer to not answer = 0",
                                        "age = 0", "Color4 = 0"), 
                                white.adjust = "hc1")

f_null <- bind_rows(existence_null, causes_null, impact_null, index_null) %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(Model = c("Existence", "Causes", "Impacts", "Belief Index")) %>%
  dplyr::select(Model, "Res.Df", "Df", "F", "Pr(>F)")

kable(f_null, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",6)))

#------------------------------------------------------------------------------#
# 3. Benjamini-Hochberg Adjusted p-values  ############################
#------------------------------------------------------------------------------#

existence_h <- feols(bf ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data = Base, weights = ~ ponde2, vcov = "hetero")

existence <- feols(bf ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                     Female + Escolaridade_3 + Religion + Income + age + Color4, 
                   data = Base, weights = ~ ponde2)

causes_h <- feols(cau ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                    Female + Escolaridade_3 + Religion + Income + age + Color4, 
                  data = Base, weights = ~ ponde2, vcov = "hetero")

causes <- feols(cau ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                  Female + Escolaridade_3 + Religion + Income + age + Color4, 
                data = Base, weights = ~ ponde2)

impacts_h <- feols(imp ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                     Female + Escolaridade_3 + Religion + Income + age + Color4, 
                   data = Base, weights = ~ ponde2, vcov = "hetero")

impacts <- feols(imp ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                   Female + Escolaridade_3 + Religion + Income + age + Color4, 
                 data = Base, weights = ~ ponde2)

index_h <- feols(index_std ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                   Female + Escolaridade_3 + Religion + Income + age + Color4, 
                 data = Base, weights = ~ ponde2, vcov = "hetero")

index <- feols(index_std ~ sk + ok + sc + ts + nep + ii + ei + pe + pi_lf + pi_cs +
                 Female + Escolaridade_3 + Religion + Income + age + Color4, 
               data = Base, weights = ~ ponde2)

# Existence
existence_std <- as.data.frame(existence_h$coeftable) %>%
  select('Pr(>|t|)') %>%
  mutate_if(is.numeric, round, digits = 5)

existence_bh <- as.data.frame(p.adjust(existence$coeftable$`Pr(>|t|)`,
                                       method = "BH")) %>%
  mutate_if(is.numeric, round, digits = 4)

existence_p <- existence_std %>%
  bind_cols(existence_bh) 

names(existence_p)[1] = "standard"
names(existence_p)[2] = "adjusted"

existence_p <- existence_p %>%
  mutate(signif = case_when(
    standard <= 0.01 ~ "***",
    standard > 0.01 & standard <= 0.05 ~ "**",
    standard > 0.05 & standard <= 0.1 ~ "*",
    standard > 0.1 ~ "-"),
    signif_adj = case_when(
      adjusted <= 0.01 ~ "***",
      adjusted > 0.01 & adjusted <= 0.05 ~ "**",
      adjusted > 0.05 & adjusted <= 0.1 ~ "*",
      adjusted > 0.1 ~ "-")) %>%
  select(standard, signif, adjusted, signif_adj)

# Causes
causes_std <- as.data.frame(causes_h$coeftable) %>%
  select('Pr(>|t|)') %>%
  mutate_if(is.numeric, round, digits = 4)

causes_bh <- as.data.frame(p.adjust(causes$coeftable$`Pr(>|t|)`,
                                    method = "BH")) %>%
  mutate_if(is.numeric, round, digits = 4)

causes_p <- causes_std %>%
  bind_cols(causes_bh) 

names(causes_p)[1] = "standard"
names(causes_p)[2] = "adjusted"

causes_p <- causes_p %>%
  mutate(signif = case_when(
    standard <= 0.01 ~ "***",
    standard > 0.01 & standard <= 0.05 ~ "**",
    standard > 0.05 & standard <= 0.1 ~ "*",
    standard > 0.1 ~ "-"),
    signif_adj = case_when(
      adjusted <= 0.01 ~ "***",
      adjusted > 0.01 & adjusted <= 0.05 ~ "**",
      adjusted > 0.05 & adjusted <= 0.1 ~ "*",
      adjusted > 0.1 ~ "-")) %>%
  select(standard, signif, adjusted, signif_adj)

# Impacts
impacts_std <- as.data.frame(impacts_h$coeftable) %>%
  select('Pr(>|t|)') %>%
  mutate_if(is.numeric, round, digits = 4)

impacts_bh <- as.data.frame(p.adjust(impacts$coeftable$`Pr(>|t|)`,
                                     method = "BH")) %>%
  mutate_if(is.numeric, round, digits = 4)

impacts_p <- impacts_std %>%
  bind_cols(impacts_bh) 

names(impacts_p)[1] = "standard"
names(impacts_p)[2] = "adjusted"

impacts_p <- impacts_p %>%
  mutate(signif = case_when(
    standard <= 0.01 ~ "***",
    standard > 0.01 & standard <= 0.05 ~ "**",
    standard > 0.05 & standard <= 0.1 ~ "*",
    standard > 0.1 ~ "-"),
    signif_adj = case_when(
      adjusted <= 0.01 ~ "***",
      adjusted > 0.01 & adjusted <= 0.05 ~ "**",
      adjusted > 0.05 & adjusted <= 0.1 ~ "*",
      adjusted > 0.1 ~ "-")) %>%
  select(standard, signif, adjusted, signif_adj)

# Index
index_std <- as.data.frame(index_h$coeftable) %>%
  select('Pr(>|t|)') %>%
  mutate_if(is.numeric, round, digits = 4)

index_bh <- as.data.frame(p.adjust(index$coeftable$`Pr(>|t|)`,
                                   method = "BH")) %>%
  mutate_if(is.numeric, round, digits = 4)

index_p <- index_std %>%
  bind_cols(index_bh) 

names(index_p)[1] = "standard"
names(index_p)[2] = "adjusted"

index_p <- index_p %>%
  mutate(signif = case_when(
    standard <= 0.01 ~ "***",
    standard > 0.01 & standard <= 0.05 ~ "**",
    standard > 0.05 & standard <= 0.1 ~ "*",
    standard > 0.1 ~ "-"),
    signif_adj = case_when(
      adjusted <= 0.01 ~ "***",
      adjusted > 0.01 & adjusted <= 0.05 ~ "**",
      adjusted > 0.05 & adjusted <= 0.1 ~ "*",
      adjusted > 0.1 ~ "-")) %>%
  select(standard, signif, adjusted, signif_adj)

# Creating table

adjusted_p <- bind_cols(existence_p, causes_p, impacts_p, index_p)

kable(adjusted_p, "latex",  booktabs = T, digits = 4, align = c(rep("c",8)))

#------------------------------------------------------------------------------#
# 4. Hypothesis Testing: Linear combination  ############################
#------------------------------------------------------------------------------#

## Existence ####

### Significant x Significant variables ####

# Subjective x Objective Knowledge
linearHypothesis(existence, c("sk - ok = 0"), 
                 white.adjust = "hc1") 

# Subjective x Scientific Consensus
sc1 <- linearHypothesis(existence, c("sk - sc = 0"), 
                 white.adjust = "hc1") 

# Subjective x NEP
linearHypothesis(existence, c("sk - nep = 0"), 
                 white.adjust = "hc1") 

# Subjective x Egalitarianism worldview
linearHypothesis(existence, c("sk - ei = 0"), 
                 white.adjust = "hc1") 

# Subjective x Personal experience
linearHypothesis(existence, c("sk - pe = 0"), 
                 white.adjust = "hc1") 

# Objective x Scientific Consensus
sc2 <- linearHypothesis(existence, c("ok - sc = 0"), 
                 white.adjust = "hc1") 

# Objective x NEP
linearHypothesis(existence, c("ok - nep = 0"), 
                 white.adjust = "hc1") 

# Objective x Egalitarianism worldview
linearHypothesis(existence, c("ok - ei = 0"), 
                 white.adjust = "hc1") 

# Objective x Personal Experience
linearHypothesis(existence, c("ok - pe = 0"), 
                 white.adjust = "hc1") 

# Scientific Consensus x NEP
sc3 <- linearHypothesis(existence, c("sc - nep = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x EI
sc4 <- linearHypothesis(existence, c("sc - ei = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x PE
sc5 <- linearHypothesis(existence, c("sc - pe = 0"), 
                                white.adjust = "hc1") 

# NEP x Personal Ei
linearHypothesis(existence, c("nep - ei = 0"), 
                                    white.adjust = "hc1")

# NEP x Personal Experience
linearHypothesis(existence, c("nep - pe = 0"), 
                                    white.adjust = "hc1")

# EI x Personal Experience
linearHypothesis(existence, c("ei - pe = 0"), 
                                    white.adjust = "hc1")

sc <-  bind_rows(sc1, sc2, sc3, sc4, sc5) %>%
  drop_na() %>%
  mutate(Variable = c("Subjective Knowledge", "Objective Knowledge", "The New Ecological Paradigm (NEP)", 
                      "Egalitarianism worldview", "Personal experience (extreme weather events"),
         .before = Res.Df)

kable(sc, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",4)))

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(existence, vcov = vcovHC(existence, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(existence, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Age (Years)" = "age",
                                                  "Trust in scientists" = "ts",
                                                  "Education: High school or equivalent" = "Escolaridade_3High school or equivalent", 
                                                  "Education: Undergraduate or more"  = "Escolaridade_3Undergraduate or more",
                                                  "Female" = "FemaleYes",
                                                  "Individualism worldview" = "ii",
                                                  "Income: 10 minimum wages or more" = "Income10 minimum wages or more",
                                                  "Income: 2 - 3 minimum wages" = "Income2 - 3 minimum wages",
                                                  "Income: 3 - 5 minimum wages" = "Income3 - 5 minimum wages",
                                                  "Income: 5 - 10 minimum wages" = "Income5 - 10 minimum wages",
                                                  "Income: Do not know/ Prefer to not answer" = "IncomeDo not know/ Prefer to not answer",
                                                  "Political ideology: Left" = "pi_lf",
                                                  "Religion: Catholic" = "ReligionCatholic",
                                                  "Religion: Evangelical Traditional" = "ReligionEvangelical Traditional",
                                                  "Religion: Others/No Relig." = "ReligionOthers/No Relig.",
                                                  "Religion: Evangelical Pentecostal or other evangelical" = "ReligionEvangelical Pentecostal or other evangelical")


levels(test_results$`Signif. Variable`) <- list("Const" = "(Intercept)",
                                                "Race: Black" = "Color4",
                                                "Egalitarianism worldview" = "ei",
                                                "Income: 1 - 2 minimum wages" = "Income1 - 2 minimum wages",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Objective knowledge" = "ok",
                                                "Personal experience (extreme weather events)" = "pe",
                                                "Political ideology: Progressive" = "pi_cs",
                                                "Scientific consensus" = "sc",
                                                "Subjective knowledge" = "sk")

# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" ))

(31+66)/(52+11+31+66) # 60.63%

# Applying Bonferroni adjustment
test_results <- test_results %>%
  mutate(adf_pvalues_bh = p.adjust(`Pr(>F)`, method = "BH"))

# Count number of significant differences at the 95% level after Bonferroni correction
test_results <- test_results %>%
  mutate(Check_II = ifelse(adf_pvalues_bh <= 0.01, "***",
                           ifelse(adf_pvalues_bh > 0.01 & adf_pvalues_bh <= 0.05, "**",
                                  ifelse(adf_pvalues_bh > 0.05 & adf_pvalues_bh <= 0.10, "*", ""))))

## Causes ####

### Significant x Significant variables ####

# Subjective x Objective Knowledge
linearHypothesis(causes, c("sk - ok = 0"), 
                                white.adjust = "hc1") 

# Objective x Scientific Consensus
linearHypothesis(causes, c("ok - sc = 0"), 
                                white.adjust = "hc1") 

# Objective x Trust in Scientists
linearHypothesis(causes, c("ok - ts = 0"), 
                                white.adjust = "hc1") 

# Objective x NEP
linearHypothesis(causes, c("ok - nep = 0"), 
                                white.adjust = "hc1") 

# Objective x Individ
ind1 <- linearHypothesis(causes, c("ok + ii = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x Trust in Scientists
linearHypothesis(causes, c("sc - ts = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x NEP
linearHypothesis(causes, c("sc - nep = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x Individ
ind2 <- linearHypothesis(causes, c("sc + ii = 0"), 
                                white.adjust = "hc1") 

# Trust in Scientists x NEP
linearHypothesis(causes, c("ts - nep = 0"), 
                                white.adjust = "hc1") 

# Trust in Scientists x II
ind3 <- linearHypothesis(causes, c("ts + ii = 0"), 
                                white.adjust = "hc1") 

# NEP x II
ind4 <- linearHypothesis(causes, c("nep + ii = 0"), 
                         white.adjust = "hc1") 

ind <-  bind_rows(ind1, ind2, ind3, ind4) %>%
  drop_na() %>%
  mutate(Variable = c("Objective Knowledge", "Scientific Consensus", "Trust in Scientists", 
                      "The New Ecological Paradigm (NEP)"),
         .before = Res.Df)

kable(ind, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",4)))

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(causes, vcov = vcovHC(causes, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(causes, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Age (Years)" = "age",
                                                  "Race: Black" = "Color4",
                                                  "Egalitarianism worldview" = "ei",
                                                  "Female" = "FemaleYes",
                                                  "Subjective knowledge" = "sk",
                                                  "Income: 1 - 2 minimum wages" = "Income1 - 2 minimum wages", 
                                                  "Income: 2 - 3 minimum wages" = "Income2 - 3 minimum wages",
                                                  "Income: 3 - 5 minimum wages" = "Income3 - 5 minimum wages",
                                                  "Income: 5 - 10 minimum wages" = "Income5 - 10 minimum wages",
                                                  "Income: 10 minimum wages or more" = "Income10 minimum wages or more",
                                                  "Income: Do not know/ Prefer to not answer" = "IncomeDo not know/ Prefer to not answer",
                                                  "Personal experience (extreme weather events)" = "pe",
                                                  "Political ideology: Left" = "pi_lf",
                                                  "Religion: Catholic" = "ReligionCatholic",
                                                  "Religion: Evangelical Traditional" = "ReligionEvangelical Traditional",
                                                  "Religion: Evangelical Pentecostal or other evangelical" = "ReligionEvangelical Pentecostal or other evangelical",
                                                  "Religion: Others/No Relig." = "ReligionOthers/No Relig.")

levels(test_results$`Signif. Variable`) <- list("Education: High school or equivalent" = "Escolaridade_3High school or equivalent", 
                                                "Education: Undergraduate or more"  = "Escolaridade_3Undergraduate or more",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Individualism worldview" = "ii",
                                                "Objective knowledge" = "ok",
                                                "Political ideology: Progressive" = "pi_cs",
                                                "Scientific consensus" = "sc",
                                                "Trust in scientists" = "ts",
                                                "Const" = "(Intercept)")
 

# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" ))

(41+52)/(42+18+41+52) # 60.78%

## Consequences ####

### Significant x Significant variables ####

# Subjective x Objective Knowledge
linearHypothesis(impact, c("sk - ok = 0"), 
                                white.adjust = "hc1") 

# Objective Knowledge x Scientific Consensus
linearHypothesis(impact, c("ok - sc = 0"), 
                                white.adjust = "hc1") 

# Objective Knowledge x NEP
linearHypothesis(impact, c("ok - nep = 0"), 
                                white.adjust = "hc1") 

# Objective Knowledge x Individualism Worldview
ii1 <- linearHypothesis(impact, c("ok + ii = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x NEP
linearHypothesis(impact, c("sc - nep = 0"), 
                                white.adjust = "hc1") 

# Scientific Consensus x Individualism Worldview
ii2 <- linearHypothesis(impact, c("sc + ii = 0"), 
                                white.adjust = "hc1") 

# NEP x Individualism Worldview
ii3 <- linearHypothesis(impact, c("nep + ii = 0"), 
                                white.adjust = "hc1") 

ii <-  bind_rows(ii1, ii2, ii3) %>%
  drop_na() %>%
  mutate(Variable = c("Objective Knowledge", "Scientific Consensus", 
                      "The New Ecological Paradigm (NEP)"),
         .before = Res.Df)

kable(ii, "latex",  booktabs = T, digits = 2, align = c("l", rep("c",4)))

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(impact, vcov = vcovHC(impact, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(impact, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Political ideology: Progressive" = "pi_cs",
                                                  "Subjective knowledge" = "sk", 
                                                  "Trust in scientists" = "ts",
                                                  "Egalitarianism worldview" = "ei",
                                                  "Personal experience (extreme weather events)" = "pe",
                                                  "Female" = "FemaleYes",
                                                  "Religion: Catholic" = "ReligionCatholic",
                                                  "Religion: Evangelical Traditional" = "ReligionEvangelical Traditional",
                                                  "Religion: Others/No Relig." = "ReligionOthers/No Relig.", 
                                                  "Income: 5 - 10 minimum wages" = "Income5 - 10 minimum wages",
                                                  "Income: Do not know/ Prefer to not answer" = "IncomeDo not know/ Prefer to not answer",
                                                  "Race: Black" = "Color4",
                                                  "Age (Years)" = "age",
                                                  "Const" = "(Intercept)")

levels(test_results$`Signif. Variable`) <- list("Political ideology: Left" = "pi_lf",
                                                "Objective knowledge" = "ok",  
                                                "Scientific consensus" = "sc",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Individualism worldview" = "ii",
                                                "Education: High school or equivalent" = "Escolaridade_3High school or equivalent", 
                                                "Education: Undergraduate or more"  = "Escolaridade_3Undergraduate or more",
                                                "Religion: Evangelical Pentecostal or other evangelical" = "ReligionEvangelical Pentecostal or other evangelical",
                                                "Income: 1 - 2 minimum wages" = "Income1 - 2 minimum wages", 
                                                "Income: 2 - 3 minimum wages" = "Income2 - 3 minimum wages",
                                                "Income: 3 - 5 minimum wages" = "Income3 - 5 minimum wages",
                                                "Income: 10 minimum wages or more" = "Income10 minimum wages or more")


# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" ))

(23+77)/(59+9+23+77) # 59.52%

#------------------------------------------------------------------------------#
# 5 Hypothesis Testing: Linear combination with standardized variables #######
#------------------------------------------------------------------------------#

# Standardizing independent variables
Base_std <- Base %>%
  fastDummies::dummy_cols(select_columns = c("Female", "Escolaridade_3", "Religion", "Income", "Color4"), remove_selected_columns = TRUE) %>%
  mutate(across(all_of(195:223), scale)) %>%
  select(-c("Female_No", "Escolaridade_3_Elementary (Primary) or less", "Religion_Atheist",
            "Income_0 - 1 minimum wages", "Color4_0"))
  
#Running the models
independent_vars <- Base_std[, 195:218]

# Existence
dataset_subset <- cbind(dependent_variable = Base_std$bf, independent_vars)

existence_std <- lm(dependent_variable ~ . + Base$age,
            data = dataset_subset,
            weights = Base_std$ponde2)

summary(existence_std)

# Causes
dataset_subset <- cbind(dependent_variable = Base_std$cau, independent_vars)

causes_std <- lm(dependent_variable ~ . + Base$age,
                    data = dataset_subset,
                    weights = Base_std$ponde2)

# Impacts
dataset_subset <- cbind(dependent_variable = Base_std$imp, independent_vars)

impact_std <- lm(dependent_variable ~ . + Base$age,
                    data = dataset_subset,
                    weights = Base_std$ponde2)

## Existence ####

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(existence_std, vcov = vcovHC(existence_std, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(existence_std, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Age (Years)" = "Base$age",
                                                  "Trust in scientists" = "ts",
                                                  "Education: High school or equivalent" = "`Escolaridade_3_High school or equivalent`", 
                                                  "Education: Undergraduate or more"  = "`Escolaridade_3_Undergraduate or more`",
                                                  "Female" = "Female_Yes",
                                                  "Individualism worldview" = "ii",
                                                  "Income: 10 minimum wages or more" = "`Income_10 minimum wages or more`",
                                                  "Income: 2 - 3 minimum wages" = "`Income_2 - 3 minimum wages`",
                                                  "Income: 3 - 5 minimum wages" = "`Income_3 - 5 minimum wages`",
                                                  "Income: 5 - 10 minimum wages" = "`Income_5 - 10 minimum wages`",
                                                  "Income: Do not know/ Prefer to not answer" = "`Income_Do not know/ Prefer to not answer`",
                                                  "Political ideology: Left" = "pi_lf",
                                                  "Religion: Catholic" = "Religion_Catholic",
                                                  "Religion: Evangelical Traditional" = "`Religion_Evangelical Traditional`",
                                                  "Religion: Others/No Relig." = "`Religion_Others/No Relig.`",
                                                  "Religion: Evangelical Pentecostal or other evangelical" = "`Religion_Evangelical Pentecostal or other evangelical`")


levels(test_results$`Signif. Variable`) <- list("Const" = "(Intercept)",
                                                "Race: Black" = "Color4_1",
                                                "Egalitarianism worldview" = "ei",
                                                "Income: 1 - 2 minimum wages" = "`Income_1 - 2 minimum wages`",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Objective knowledge" = "ok",
                                                "Personal experience (extreme weather events)" = "pe",
                                                "Political ideology: Progressive" = "pi_cs",
                                                "Scientific consensus" = "sc",
                                                "Subjective knowledge" = "sk")

# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" ))

# (83+27) / (14+36+27+83) = 69%

# Applying Bonferroni adjustment
test_results <- test_results %>%
  mutate(adf_pvalues_bh = p.adjust(`Pr(>F)`, method = "BH"))

# Count number of significant differences at the 95% level after bonferroni correction
test_results <- test_results %>%
  mutate(Check_II = ifelse(adf_pvalues_bh <= 0.01, "***",
                           ifelse(adf_pvalues_bh > 0.01 & adf_pvalues_bh <= 0.05, "**",
                                  ifelse(adf_pvalues_bh > 0.05 & adf_pvalues_bh <= 0.10, "*", ""))))

## Causes ####

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(causes_std, vcov = vcovHC(causes_std, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(causes_std, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Age (Years)" = "Base$age",
                                                  "Race: Black" = "Color4_1",
                                                  "Egalitarianism worldview" = "ei",
                                                  "Female" = "Female_Yes",
                                                  "Subjective knowledge" = "sk",
                                                  "Income: 1 - 2 minimum wages" = "`Income_1 - 2 minimum wages`", 
                                                  "Income: 2 - 3 minimum wages" = "`Income_2 - 3 minimum wages`",
                                                  "Income: 3 - 5 minimum wages" = "`Income_3 - 5 minimum wages`",
                                                  "Income: 5 - 10 minimum wages" = "`Income_5 - 10 minimum wages`",
                                                  "Income: 10 minimum wages or more" = "`Income_10 minimum wages or more`",
                                                  "Income: Do not know/ Prefer to not answer" = "`Income_Do not know/ Prefer to not answer`",
                                                  "Personal experience (extreme weather events)" = "pe",
                                                  "Political ideology: Left" = "pi_lf",
                                                  "Religion: Catholic" = "Religion_Catholic",
                                                  "Religion: Evangelical Traditional" = "`Religion_Evangelical Traditional`",
                                                  "Religion: Evangelical Pentecostal or other evangelical" = "`Religion_Evangelical Pentecostal or other evangelical`",
                                                  "Religion: Others/No Relig." = "`Religion_Others/No Relig.`")

levels(test_results$`Signif. Variable`) <- list("Education: High school or equivalent" = "`Escolaridade_3_High school or equivalent`", 
                                                "Education: Undergraduate or more"  = "`Escolaridade_3_Undergraduate or more`",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Individualism worldview" = "ii",
                                                "Objective knowledge" = "ok",
                                                "Political ideology: Progressive" = "pi_cs",
                                                "Scientific consensus" = "sc",
                                                "Trust in scientists" = "ts",
                                                "Const" = "(Intercept)")


# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" )) 

# (78 + 32) / (27 + 16 + 78 +32) = 72%

## Consequences ####

### F-Test: Significant x Insignificant variables ####

# Extract the p-values of the coefficients from the summary
p_values <- coeftest(impact_std, vcov = vcovHC(impact_std, type = "HC1"))[, "Pr(>|t|)"]

# Create a vector of the significant variables
significant_vars <- names(p_values[p_values <= 0.05])

# Create a vector of the insignificant variables
insignificant_vars <- names(p_values[p_values > 0.05])

# Use the "linearHypothesis()" function to test for linear combinations of the coefficients
test_results <- data.frame() # Create an empty data frame to store the results
for (i in 1:length(significant_vars)) {
  for (j in 1:length(insignificant_vars)) {
    test <- linearHypothesis(impact_std, paste(significant_vars[i], "-", insignificant_vars[j], "=0"), white.adjust = "hc1")
    test_results <- rbind(test_results, c(insignificant_vars[j], significant_vars[i], test$`Res.Df`[2], test$`D`[2], test$F[2], test$`Pr(>F)`[2]))
  }
}

# Name the columns of the results data frame
colnames(test_results) <- c("Insignif. Variable", "Signif. Variable", "Res.Df", "Df", "F", "Pr(>F)")

# Converting columns to numeric
test_results <- test_results %>%
  mutate_at(vars(`Res.Df`, `Df`, `F`, `Pr(>F)`), as.numeric) %>%
  mutate(Check = ifelse(`Pr(>F)` <= 0.01, "***",
                        ifelse(`Pr(>F)` > 0.01 & `Pr(>F)` <= 0.05, "**",
                               ifelse(`Pr(>F)` > 0.05 & `Pr(>F)` <= 0.10, "*", ""))),
         `Insignif. Variable` = as.factor(`Insignif. Variable`),
         `Signif. Variable` = as.factor(`Signif. Variable`)) 

levels(test_results$`Insignif. Variable`) <- list("Political ideology: Progressive" = "pi_cs",
                                                  "Subjective knowledge" = "sk", 
                                                  "Trust in scientists" = "ts",
                                                  "Egalitarianism worldview" = "ei",
                                                  "Personal experience (extreme weather events)" = "pe",
                                                  "Female" = "Female_Yes",
                                                  "Religion: Catholic" = "Religion_Catholic",
                                                  "Religion: Evangelical Traditional" = "`Religion_Evangelical Traditional`",
                                                  "Religion: Others/No Relig." = "`Religion_Others/No Relig.`", 
                                                  "Income: 5 - 10 minimum wages" = "`Income_5 - 10 minimum wages`",
                                                  "Income: Do not know/ Prefer to not answer" = "`Income_Do not know/ Prefer to not answer`",
                                                  "Race: Black" = "Color4_1",
                                                  "Age (Years)" = "Base$age",
                                                  "Const" = "(Intercept)")

levels(test_results$`Signif. Variable`) <- list("Political ideology: Left" = "pi_lf",
                                                "Objective knowledge" = "ok",  
                                                "Scientific consensus" = "sc",
                                                "The New Ecological Paradigm (NEP)" = "nep",
                                                "Individualism worldview" = "ii",
                                                "Education: High school or equivalent" = "`Escolaridade_3_High school or equivalent`", 
                                                "Education: Undergraduate or more"  = "`Escolaridade_3_Undergraduate or more`",
                                                "Religion: Evangelical Pentecostal or other evangelical" = "`Religion_Evangelical Pentecostal or other evangelical`",
                                                "Income: 1 - 2 minimum wages" = "`Income_1 - 2 minimum wages`", 
                                                "Income: 2 - 3 minimum wages" = "`Income_2 - 3 minimum wages`",
                                                "Income: 3 - 5 minimum wages" = "`Income_3 - 5 minimum wages`",
                                                "Income: 10 minimum wages or more" = "`Income_10 minimum wages or more`")


# Print the results as a LaTeX table
kable(test_results, "latex",  booktabs = T, digits = 2, align = c(rep("l",2), rep("c",5)))

# Count number of significant differences at the 95% level
test_results %>%
  group_by(Check) %>%
  summarise(n_sign = sum(Check == "***" | Check == "**" ),
            n_insig = sum(Check == "*" | Check == "" ))

# (93+23)/(93+23+46+7) = 69%

#------------------------------------------------------------------------------#
# 6. Step wise: Only significant coefficients + add other variables ######
#------------------------------------------------------------------------------#

## Existence #### 

step_existence <- feols(bf ~ csw(sk, ok, sc, nep, ei, pe, pi_cs, Color4, #Significant
                                 ts, ii, pi_lf, Female, Escolaridade_3, Religion, Income, age), #Insignificant 
                        data = Base, 
                        weights = ~ ponde2, 
                        vcov = "hetero")

names(step_existence) <- c(1:16)

step_existence <- step_existence[c(8,9,10,11,12,13,14,15,16)]

names(step_existence) <- c(1:9)

etable(step_existence,
       title = "OLS Results - Existence of Climate Change",
       digits = 3,
       dict = mydict,
       signif.code = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/stepwise_existence.tex")

modelsummary(step_existence,
             coef_omit = 'Interc',
             estimate = "{estimate}{stars} ({std.error})",
             statistic =  NULL,
             vcov = "robust",
             coef_rename = mydict,
             gof_map = c("nobs", "r.squared"),
             stars = c('*' = .10, '**' = .05, '***' = .01),
             output = "kableExtra")

## Causes ####

step_causes <- feols(cau ~ csw(ok, sc, ts, nep, ii, pi_cs, Escolaridade_3, #Significant
                               sk, ei, pe, pi_lf, Female, Religion, Income, age, Color4), #Insignificant 
                        data = Base, 
                        weights = ~ ponde2, 
                        vcov = "hetero")

names(step_causes) <- c(1:16)

step_causes <- step_causes[c(7,8,9,10,11,12,13,14,15,16)]

names(step_causes) <- c(1:10)

etable(step_causes,
       title = "OLS Results - Anthrogenic Causes of Climate Change",
       digits = 3,
       dict = mydict,
       signif.code = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/step_causes.tex")

modelsummary(step_causes,
             coef_omit = 'Interc',
             estimate = "{estimate}{stars} ({std.error})",
             statistic =  NULL,
             vcov = "robust",
             coef_rename = mydict,
             gof_map = c("nobs", "r.squared"),
             stars = c('*' = .10, '**' = .05, '***' = .01),
             output = "kableExtra")

## Consequences ####

step_imp <- feols(imp ~ csw(ok, sc, nep, ii, pi_lf, Escolaridade_3, #Significant
                            Religion, Income, sk, ts, ei, pe, pi_cs, Female, age, Color4), #Insignificant 
                     data = Base, 
                     weights = ~ ponde2, 
                     vcov = "hetero")

names(step_imp) <- c(1:16)

step_imp <- step_imp[c(6,7,8,9,10,11,12,13,14,15,16)]

names(step_imp) <- c(1:11)

etable(step_imp,
       title = "OLS Results - Consequences of Climate Change",
       digits = 3,
       dict = mydict,
       signif.code = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/step_consequences.tex")

modelsummary(step_imp,
             coef_omit = 'Interc',
             estimate = "{estimate}{stars} ({std.error})",
             statistic =  NULL,
             vcov = "robust",
             coef_rename = mydict,
             gof_map = c("nobs", "r.squared"),
             stars = c('*' = .10, '**' = .05, '***' = .01),
             output = "kableExtra")
